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1.  INTRODUCTION 


w  Many  experimental  situations  require  consideration  of  a  large 
number  of  factors.  In  such  situations  there  is  often  a  need,  because 
of  resource  limitations,  for  an  efficient  method  of  factor  screening. 

One  such  method  is  two-stage  group  screening.  In  this  method  the  in¬ 
dividual  factors  (each  at  two  levels)  are  partitioned  into  groups,  for¬ 
ming  group  factors.  By  assigning  the  same  level  to  all  component  fac¬ 
tors  within  each  group,  the  group  factors  are  tested  as  if  they  were 
single  factors.  All  factors  in  groups  found  to  have  a  significant  ef¬ 
fect  are  then  tested  individually  in  a  second-stage  experiment. 

In  Desmatics  Technical  Report  No.  113-9  (August  1983),  "On  the 
Performance  of  Two-Stage  Group  Screening  Experiments,"  a  methodology 
la  developed  with  which  to  evaluate  the  performance  of  the  two-stage 
group  screening  procedure.  Also  in  that  report,  a  computer-aided  search 
routine  is  described.  This  routine  may  be  used  to  help  select  a  satis¬ 
factory  group  screening  strategy.  In  the  present  report,  we  review  the 
basic  screening  model,  give  a  brief  description  of  the  search  routine, 
and  provide  a  listing  of  the  program.  The  search  program,  which  is 
written  in  standard  FORTRAN,  is  easy  to  use,  has  intuitive  appeal,  and 
provides  practical  insight  into  the  use  and  selection  of  a  two-stage 


group  screening  strategy 


MODEL  ASSUMPTIONS 


We  assume  the  following  model: 


*1  ■  60  +  jLVu  +  ^ 


where  K  denotes  the  number  of  factors  to  be  screened  and  =  ±1  is 
the  level  of  the  j**1  factor  in  the  i**1  run.  We  additionally  assume: 

1.  k  _>  1  (k  unknown)  of  the  K  factors  are  active  (i.e.,  have 
a  nonzero  effect)  and  the  remaining  (K-k)  are  inactive, 

2.  all  k  active  factors  have  the  same  absolute  effect,  A  >  0, 


3.  the  error  terms  {e^}  are  i.i.d.  N(0,  a2)  random  variables 
with  a2  unknown. 

In  view  of  assumptions  (1)  and  (2)  we  let  ji(i)  for  i  ■  0,  1,  2,  ... 
denote  the  case  in  which  i  active  effects  are  equal  to  -A,  (k-i)  active 

effects  are  equal  to  +A,  and  the  (K-k)  inactive  effects  are  equal  to 
zero.  The  3(0)  and  jJ(k)  cases  correspond  to  the  situation  in  which  all 
active  effects  are  in  the  same  direction.  In  these  cases  there  can  be 
no  cancellation  of  effects  within  a  group.  For  j3(i)  cases  with  i  i*  0 
and  i  +  k  cancellation  can  occur.  The  search  routine  which  we  have  de¬ 
veloped  considers  the  3^(0)  and  ji([k/2])  cases,  the  latter  being  the 
case  in  which  the  chance  of  cancellation  is  greatest,  among  all  ji(i) 
cases. 


In  the  two-stage  group  screening  strategy  considered  here,  we 


assume  that  the  K  factors  are  partitioned  randomly  into  G  groups  of 
size  g;  if  K  is  not  a  multiple  of  g,  we  assume  that  the  group  sizes 
are  taken  as  "evenly"  as  possible.  To  analyze  the  results  of  the  first 


and  second  stages,  we  use  the  multifactorial  designs  of  Plackett  and 
Burman  (PB).  We  denote  the  levels  of  the  significance  tests  performed 
at  the  end  of  the  first  and  second  stages  by  and  o^,  respectively. 
Following  the  notation  used  in  Technical  Report  No.  113-9,  we  denote 
such  a  strategy  by  GS(  g,  a.,  a»). 


3.  PROGRAM  DESCRIPTION 


We  can  define  three  separate  measures  of  performance.  These 


are: 


E 


U 


Power .  We  denote  by  A  the  number  of  active  factors  that  are 
detected  correctly,  and  we  define  EA  «  100E(A)/k  as  a  per¬ 
centage  measure  of  the  power  of  a  GS  strategy  to  detect  active 
factors. 

Type  I  Error.  We  denote  by  U  the  number  of  inactive  factors 
that  are  declared  active,  and  we  define  Ey  -  100E(U)/(K-k)  as 
a  percentage  measure  of  Type  1  error. 

Relative  Testing  Cost.  We  denote  by  R  the  total  number  of  runs 
required  by  a  GS  strategy,  and  we  define  ER  -  100E(R) /B(K+1)  as 
a  percentage  measure  of  expected  relative  testing  cost  where 
B (K+l )  denotes  the  number  of  runs  required  by  a  PB  design  for 
(K+l)  factors.1 

In  general,  an  increase  in  or  will  increase  both  and 


An  increase  in  ct^  will  also  Increase  ER,  although 


does  not 


depend  on  a^.  It  is  clear  that  trade-offs  between  EA,  Ey,  and  Ej^  must 
be  made  by  the  experimenter  in  order  to  select  a  "good"  GS(g,  a^,  a^) 
strategy.  The  appropriate  trade-offs  would  most  likely  be  determined 
by  economic  and  operational  considerations. 

_ A  general  practical  approach  to  this  problem  is  to  maximize 

^B  designs  exist  for  numbers  of  runs  that  are  multiples  of  four 
only.  Accordingly,  one  can  define  B(x)  ■  x+4-x(mod  4),  We  use  x  “  K+l 
instead  of  x  ■  K  in  our  definition  of  E^  to  guarantee  at  least  one 
degree  of  freedom  to  estimate  error. 
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power  (E^)  subject  to  fixed  constraints  on  relative  testing  cost  (E^) 
and  Type  I  Error  (Ey).  This  approach  is  the  basis  of  the  search  routine. 
That  is,  in  order  to  apply  the  search  program,  the  user  must  first  input 
values  of  K,  k,  and  A/a,  the  signal-to-noise  ratio.  In  addition,  the 
user  must  specify  a  maximum  tolerable  relative  testing  cost,  say  E*  , 

E*  and 

Ey  _<  E*,  the  search  program  determines,  for  various  group  sizes,  the 
values  of  and  ct^  that  maximize  power  (E^) .  From  the  program  output, 
the  group  size  (and  corresponding  and  a ^  levels)  which  gives  the 
greatest  power  may  then  be  selected.  The  basic  steps  comprising  the 
search  algorithm  are  outlined  in  Table  1. 

The  logic  of  the  search  is  based  on  the  premise  that  power  is  an 
increasing  function  of  testing  cost.  That  is,  the  more  runs  that  are 
invested,  the  greater  is  the  power.  Consequently,  we  can  first  determine 
independently  of  a ^  since  E^  is  a  function  of  alone.  Moreover, 

■  E*.  Once  has  been  deter¬ 
mined,  it  is  not  difficult  to  determine  the  optimal  value  of  ot ^  subject 
to  Ey  <  E*.  We  remark  that  a  ’’solution"  in  which  ER  is  strictly  less 
than  E*  would  be  counterintuitive. 


we  find  (if  possible)  the  otj  such  that  E^ 


and  a  maximum  tolerable  Type  I  error,  say  E*.  Subject  to  ER  £ 


Input  values  for  K,  k,  and  A/a  . 

Input  maximum  tolerable  values  for  ER  and  Ey. 

Assume  £(0)  case  and  g»2. 

Determine  so  that  ER  attains  maximum  allowable  value 

For  the  determined  in  Step  4,  determine  the  a. 
that  maximizes  E„  subject  to  constraint  specified 
in  Step  2. 

Calculate  E^  and  Ey  for  given  GSCgta^a^)  strategy. 

Repeat  Steps  4,  5,  and  6  as  long  as  g  _<  min(8,  K/2) . 

Reset  g-2  and  repeat  Steps  4  through  7  for  g(Tk/2]) 
case. 


Table  1.  Outline  of  Two-Stage  Group  Screening 
Search  Algorithm 


4.  USER’S  GUIDE 


In  this  section  we  discuss  several  key  aspects  of  the  search  pro¬ 
gram,  detail  the  input  considerations  and  requirements,  describe  the 
various  program  modules,  and  present  some  sample  printout.  In  this  sec¬ 
tion,  numbers  appearing  in  brackets  refer  to  line  numbers  in  the  program 
listing  which  is  provided  in  Section  5. 

4.1  IMSL  Subprograms 

The  search  program  makes  direct  use  of  three  IMSL  (The  International 
Mathematical  and  Statistical  Library)  subprograms.  These  are  MDTN, 

MDSTI,  and  MDBIN.  The  IMSL  subroutine  MDTN  calculates  a  non-central 
t  probability  and  MDBIN  calculates  a  binomial  probability.  The  IMSL 
routine  MDSTI  is  used  to  obtain  inverse  values  of  Student's  t  probability 
distribution  function. 

Because  of  copyright  regulations,  however,  the  source  listings  for 
each  of  these  routines  (and  any  other  IMSL  routines  called  by  these 
routines)  are  not  included  with  the  search  program  listing  given  in 
Section  5.  Therefore,  to  apply  the  search  program  the  user  must  either 
have  calling  access  to  the  ISML  package  of  subprograms  or  provide  equi¬ 
valent  alternative  routines.  The  Appendix  contains  additional  documen¬ 
tation  of  the  MDTN,  MDSTI,  and  MDBIN  routines. 

4.2  ERRSET  Subroutine 

The  ERRSET  subroutine  is  part  of  the  Extended  Error  Handling 
Facili ty  provide<  by  IBM.  As  called  by  the  search  program,  see  [22], 
ERRSET  own  re.  as  all  underflow  messages  and  allows  for  an  unlimited 


number  of  such  occurrences.  The  standard  corrective  action  for  under¬ 
flow,  i.e.,  set  result  register  to  zero,  is  unchanged.  Moreover, 
no  traceback  is  printed.  Underflow  can  arise  in  the  search  routine 
when  tail  probabilities  that  are  essentially  zero  must  be  calculated. 

4.3  Input  Considerations 

Application  of  the  two-stage  group  screening  search  routine 
requires  a  minimal  amount  of  effort  on  the  part  of  the  user.  Only 
two  input  data  (control)  cards  must  be  prepared.  Furthermore,  to 
simplify  the  input  process  format-free  input  is  assumed. 

The  first  input  data  card  must  contain  values  for  K,  k,  and  A/a, 
in  that  order.  The  second  input  data  card  must  contain  values  for 
E*  and  E*,in  that  order.  Furthermore,  the  user  must  express  E*  and  E* 
as  percentages. 

The  values  of  K,  k,  and  A/a  are  limited  only  by  their  output 
formats,  see  [43/44].  In  addition,  for  ease  of  calculation* ,  k  is 
assumed  to  be  even.  If  k  is  odd,  the  program  automatically  redefines 
k  as  k+1  and  prints  a  message  to  that  effect. 

The  search  program  also  contains  several  checks  for  invalid 
input  parameters,  see  [27/31].  The  following  error  condition  codes 
are  used: 

Condition  Code  Description 

1  E*  >  100.  In  this  case,  we  recommend  the 

use- of  a  PB  design  (i.e.,  g*l)  for  screening 

2  E*  <  °. 

*If  k  is  even,  there  are  an  equal  number  of  positive  (+A)  and  negative 
(-A)  active  effects  in  the  ( [k/2 ) ]  case.  In  this  case  the  search  program 
can  exploit  the  existing  symmetry  in  order  to  substantially  reduce  the 
extent  of  the  required  performance  calculations.  See  Technical  Report  No.  113-9 
for  further  details. 


Condition  Code 


3 

4 

5 

4.4  Program  Modules 
A  brief  description  of  each  program  module  in  the  search  program 

is  provided  in  Table  2.  This  description  provides  information  about  the 
purpose  of  the  module,  its  location  in  the  source  listing,  and  other 
modules  which  call  and  are  called  by  it.  The  PER  subprogram  is  the 
key  subroutine  that  evaluates  performance.  PER  can  evaluate  the  per¬ 
formance  of  any  GS(g,  a^,  o^)  strategy. 

4.5  An  Example 

Table  3  contains  sample  computer  printout  from  the  search  routine 

for  the  case  K  *  60,  k  -  8,  A/a  *  2,  E*  -  50%,  and  E*  =*  5%.  Asterisks 

R  U 

appearing  in  the  printout  signify  a  group  size  in  which  the  value  of  E 

R 

at  Oj  ■  0  equals  or  exceeds  E*.  In  Table  3,  for  instance,  when  g  *  2 
the  first-stage  PB  experiment  requires  32  runs,  thus  leaving  no  runs 
available  for  the  second-stage  follow-up  experiment.  Asterisks,  there¬ 
fore  indicate  that  E*  is  out  of  range  for  the  given  group  size. 

The  search  routine  is  based  on  the  performance  results  derived  in 
Technical  Report  No.  113-9.  These  results,  however,  are  only  appli¬ 
cable  when  K  is  a  multiple  of  g.  Therefore,  for  group  sizes  where 
this  restriction  is  not  met,  the  program  redefines  K  as  the  nearest 
multiple  of  g  and  denotes  the  new  value  as  KSTAR.  Performance 
results  are  then  calculated  as  if  there  are  KSTAR  factors  to  be  screened. 
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Table  3.  Sample  Computer  Printout  From  Two-Stage  Group  Screening  Search 
Routine  When  K**60,  k«8,  and  A/o*2.  Maximum  E..  Specified  As  5% 
And  Maximum  E„  Specified  As  50%  (Equivalently,  E(R)*32  Runs). 
Note:  Power,  Type  I  Error,  and  Relative  Testing  Cost  Are 
Expressed  As  Percentages. 


For  example,  in  Table  3  for  g  ■  7,  the  number  of  groups  is  taken  as 
G  *  9,  and  thus  KSTAR  ■  gG  *  63.  It  seems  reasonable  that  the  cor¬ 
responding  results  for  g  ■  7  should  be  comparable  to  the  true  per¬ 
formance  had  the  original  K  *  60  factors  been  partitioned  as  "evenly" 
as  possible  -  i.e.,  6  groups  of  size  7  and  3  groups  of  size  6. 

An  examination  of  Table  3  shows  that  the  GS(7,  .00714,  .19609) 
strategy  has  greatest  power  in  the  j3(0)  case  and  the  GS(7»  .01312,  .18467) 
strategy  has  greatest  power  in  the  j3(4)  case.  Greater  power,  of  course, 
is  attained  in  the  J3(0)  case,  when  no  cancellation  can  occur. 

The  search  routine  presented  and  described  in  this  report  provides 
guidance  in  using  and  selecting  a  satisfactory  GS  strategy.  The  routine 
supplies  the  user  with  quantitative  Information  needed  to  determine 
whether  two-stage  group  screening  is  suitable  for  a  particular  application. 


5.  PROGRAM  LISTING 
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TWO-STAGE  GROUP  SCREENING  SEARCH  ROUTINE 


AUGUST  1983.  DESMATICS,  INC. 

P.O.  BOX  618 

STATE  COLLEGE,  PA  16801 


USER  MUST  SUPPLY  FOLLOWING  TWO  CONTROL  CARDS  (FREE  FORMAT): 
CARD  #1 

K  -  TOTAL  NUMBER  OF  FACTORS  TO  BE  SCREENED 

NACT  -  NUMBER  OF  ACTIVE (I .E  .  ,NONZ ERO )  FACTORS 
SNR  -  SIGNAL -TO -NOISE  RATIO  FOR  ACTIVE  FACTORS 

CARD  #2 

RTC  -  MAXIMUM  RELATIVE  TESTING  COST  (PER  CENT) 

TIE  -  MAXIMUM  TYPE  I  ERROR  (PER  CENT) 


PROGRAM  IS  BASED  ON  DESMATICS,  INC.  TECHNICAL  REPORT  NO. 
113-9.  ALSO  SEE  USER'S  GUIDE,  TECHNICAL  REPORT  NO.  113-10. 


DIMENSION  PR (  7 , 5  )  ,  IPR  (  7  ,2  ) 

INTEGER  G  ,GS  ,GL IM 

CALL  ERRSET  (2  08  ,2  56,-1  ,1  ,1  ,0) 

READ  VALUES  FOR  K,  NACT,  SNR 

READ  (5,*)  K  , NACT ,SNR 

READ  VALUES  FOR  RTC, TIE 

READ  (5,*)  RTC, TIE 
IF  (RTC  .GE.  100.)  STOP  1 
IF  (RTC  .LE.  0.)  STOP  2 

IF  (TIE  .LE.  0.  .OR.  TIE  .GE.  100.)  STOP  3 
IF  (NACT  .GE.  K)  STOP  4 
IF  (SNR  .LE.  0.)  STOP  5 
IPBK-K+5-M0D(K+l  ,4) 

PBK-IPBK 
ER-PBK*RTC/1 00 
1F1  -0 

IF  (M0D(NACT  ,2  )  .EQ.  0)  GO  TO  5 
IF1-1 

NACT-NACT+1 
5  CONTINUE 
WRITE  (6,10) 

10  FORMAT  ('1  ',  'TWO-STAGE  GROUP  SCREENING  SEARCH  ROUTINE') 
WRITE  (6,15)  K, NACT, SNR 

15  FORMAT  ('O', 'FOR  THE  CASE:  FACTORS  -' ,14 ,/ ,1  5X  , 'ACTIVE  -  ', 

*  13  ,  /  ,1 5X  ,'SIGNAL-TO -NOISE  RATIO  »',F4.1) 

IF  (IF1  .EQ.  0)  WRITE  (6,20) 

2  0  FORMAT  ('-') 

IF  (IF1  .EQ.  1)  WRITE  (6,25) 

2  5  FORMAT  ('O', 'NOTE:  THE  NUMBER  OF  ACTIVE  FACTORS  HAS  BEEN',/ 

*  ?X  .'INCREMENTED  BY  ONE  SO  THAT  IT  IS  EVEN') 

VRITE  (6,30)  RTC, ER, TIE 
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30  FORMAT  ('-',/ ,1 X ,  'RELATIVE  TESTING  COST  -  '.F4.1.3X, 

*  '(EXPECTED  NUMBER  OF  RUNS  -' ,F  6 . 1  ,')',/ , 

«  IX, 'MAXIMUM  TYPE  I  ERROR  -  '  ,F5.2  ,// ,1X  ,'NOTE:  ', 

*  'POWER,  TYPE  I  ERROR,  AND  RELATIVE  TESTING  COST  ARE 

*  '  EXPRESSED  AS  PERCENTAGES') 

GLIM-MIN0(K/2  ,8) 

I-NACT 
00  90  J-l  ,2 

IF  (J  .EQ.  2  )  I-NACT/2 

START  LOOP  OVER  GROUP  SIZES 
00  60  GS-2  .GLIM 
M-GS-1 

Z-FLOAT(K)/GS 

G-Z 

IF  (Z-G  .GT.  .5)  G-G+l 
KS-G*GS 

PBG-G+5-M0D(G+l  ,4) 

PBKS-KS+5-M0D(KS+l ,4) 

GAM-1  00*PBG/PBKS 

IF  ( RTC  .GT.  GAM+3  50./PBKS)  GO  TO  35 
PR(M  ,1  )-l  0 
PR(M  ,2  )-l  0 
1PR (M  ,1  )-1000 
1PR (M  ,2  )-10000 
PR  (M  ,3  )-l  000 
PR (M  ,4  )-l  00 
PR(M  ,5 )-l 00 
GO  TO  6  0 
35  CONTINUE 

DETERMINE  FIRST  STAGE  SIGNIFICANCE  LEVEL 

A-0 

B-l 

PA- GAM 
FB-1  OO+GAM 
SL1-(RTC-GAM)/1 00 

CALL  PER  (1  ,REA,REU,RER,KS,NACT,SNR,I,GS,SL1  ,SL2) 
F-RER 

40  IF  (ABS(F-RTC)  .LT.  .  02  5)  GO  TO  55 

IF  (F  .LT.  RTC)  GO  TO  45 

B-SL1 
PB-F 

GO  TO  50 
45  A-SL1 

FA-F 

5  0  CONTINUE 

SLOPE-(FB-FA)/(B-A) 

SL1 -A+(RTC-FA)/ SLOPE 

CALL  PER  (1  ,R  EA  ,  REU  , RER  ,K  S  ,N ACT  ,  SNR  ,1  ,GS  ,  SL1  ,SL2  ) 

P-RER 

GO  TO  40 


X 


O®®M»UiM*lW«O'O®Ma>uiMJM>-,ovo®N,O'^*'wNMO*0D«J»Ui»WNH,ovO»MOiWi»>wf5M 


3  5 


c 


c 


C 

C 

C 

C 


6  0 

65 

70 

75 


PR (M  ,1  )-SLl 
lPR (M  ,1  )-G 
IPR (M  ,2  )-KS 

CALL  PER  (0,REA,REU,RER,KS,NACT,SNR,I  ,GS,SL1  ,.l  ) 

DETERMINE  SECOND-STAGE  SIGNIFICANCE  LEVEL 
SL2 “AM INI ( 1  .  ,T1 E/(l 0*REU) ) 

PR  (M  ,2  )-SL2 

CALCULATE  POWER  FOR  SPECIFIED  STRATEGY 
CALL  PER  ( 0 ,REA , REU , RER  ,K  S  ,NACT , SNR , I , GS , SL1  , SL2  ) 

PR (M  ,3  )-REA 
PR (M  ,4  )-REU 
PR(M , 5 ) “RER 
CONTINUE 
I-NACT-I 
WRITE  (6,65)  I 

FORMAT  ('-',//////,  IX  ,  'SEARCH  RESULTS  FOR  B(',I2,')  CASE:') 
WRITE  (6,70) 

FORMAT  ('0' ,6X .'STRATEGY' ,1 IX , 'GROUPS ', 4X , 'KSTAR ', 4X , 

'PO  WER' , 4X , 'TYPE  I  ERROR ', 4X , 'RELATIVE  TESTING  COST') 
WRITE  (6,75) 

FORMAT  (  '  +  '  ,6X  ,  '  '.11X,'  '  ,4X  ,  '  ',4X, 

' _ ',4X/  _ '  ,4x  /_ _ ') 

DO  80  GS-2  , GLIM 
M-GS-1 


60  WRITE  (6,85)  GS  ,  (PR(M  ,  I )  ,  1-1  ,2  )  ,  ( I  PR  (M  ,  I )  ,  1-1  ,2  )  , 

-  (  PR(M , I ) , I“3  , 5 ) 

85  FORMAT  ('  '  ,  'GS  (  '  ,  II  , '  ,  '  ,F  7 . 5  ,  '  ,  '  ,F  7 . 5  ,  '  )  '  ,  5X  ,13  ,  6X  ,  1 4  , 
»  4X  , F  6  •  2  ,  7X  , F  6  •  3  , 1  6X  , F  4 . 1  ) 

90  CONTINUE 
STOP 
END 


SUBROUTINE  PER(ICODE  , REA  , REU  ,RER  ,K  ,NACT  ,R  ,11  ,GS,SL1  ,SL2  ) 

IF  ICODE  EQUALS  1  THEN  ONLY  RER  RETURNED 

INTEGER  G ,GS 
REAL  CV ( 4 ) 

REAL  *  8  C,R1  , RM  ,ES ,RA ,RB 
G-K/GS 

KNACT-K-NACT 
PBK“K+5— MOD  (K+l  ,4) 

IPB-G+5-M0D(G+l ,4) 

IDF1 ■IPB-G-1 
PB1 -IPB 
RIDF1  -IDF1 

CALL  MDSTI ( SL1  ,RIDF1  ,TVAL,IER) 

RA-R1 (K , GS ,NACT , 1 1 ,R,TVAL,IDF1  ,PBl ) 

RB-RM(K  ,GS  ,NACT  , II  ,R,TVAL,IDF1  ,PB1  ) 

ES-NACT*RA+KNACT*RB 
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151  . 

152  . 

153. 

154. 
1  55. 
156. 
1  57. 
158. 
1  59. 
160. 
161  . 
1  62  . 

163. 

164. 

165. 

166. 

167. 

168. 

1  69. 
170. 

171  . 

172  . 

173. 

174. 

175. 

176. 

177. 

178. 

179. 

180. 
181  . 
182  . 

183. 

184. 

185. 

186. 

187. 

188. 

189. 

190. 

191  . 

192  . 

193. 

194. 

195. 

196. 

197. 

198. 

199. 

2  00. 


RER-1 00. * (PB1 +3 . 5+E  S) / PBK 
1F(IC0DE.EQ.1 )  RETURN 
DA-1  . 

IF ( SL2  . EQ . 1  )  CO  TO  4 

DO  3  1-1  ,4 

Rl-I 

3  CALL  MDSTI (SL2  ,RI  ,CV(I )  ,IER) 

PG-ES/K 

DA-PI  (GS  ,G  ,PG  ,R  ,CV) 

4  CONTINUE 
KEA-1 00. *DA*RA 
KEU-1 00. *SL2 *RB 
RETURN 

END 

DOUBLE  PRECISION  FUNCTION  C(N,R) 

INTEGER  R 
C-0.  ODO 

iF(R.GT.N)  RETURN 
C-l  .  ODO 

lF(R.EQ.O.OR.R.EQ.N)  RETURN 
J-R 

IF  (R  .  GT  .  N/2  )  J-N-R 
DO  1  0  1-1  ,J 

10  C-C*  (N-I+l  )  /  (J-I+l  ) 

RETURN 

END 

DOUBLE  PRECISION  FUNCTION  R1 (K  ,GS ,NACT ,1 ,R  ,TVAL ,IDF1  , PB1 ) 
INTEGER  GS.GSJ.GSJL 
REAL *8  C  ,  PJL 
KNACT-K-NACT 
11  -1-1 

NACTI-NACT-I 
D-SQRT( PB1 )*R 
Ll  -MINO(GS  ,1) 

Rl-O.DO 
DO  1  0  J-l  , Ll 
Jl-J-1 
GSJ-GS-J 

L2  -MIN0(GSJ  , NACT I )+l 
DO  2  0  L-l  , L2 
GSJL-GSJ-L+1 

PJL-C ( 1 1  ,J1 )*C(NACTI  ,L-1 )*C(KNACT ,GSJL) 

DP-D*(L-J-1  ) 

CALL  MDTN(TVAL  , IDF1  ,DP,A1  ,IER) 

CALL  MDTN(TVAL  ,IDF1  ,-DP,A2  ,IER) 

2  0  R1-R1+PJL*(2  .D0-A1-A2  ) 

10  CONTINUE 

K1-R1/C(K-1  ,GS-1 ) 
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c 


RETURN 

END 


2  01  . 
2  02  . 
2  03. 
2  04. 
2  05. 
2  06. 
2  07. 
2  08. 
2  09. 
2  10. 
211  . 
2  12  . 
2  13. 
2  14. 

215. 

216. 
2  17. 
218. 

219. 

220. 
221. 
2  22  . 

223. 

224. 

225. 

226. 
227. 
22  8. 
22  9. 
230. 
2  31  . 
2  32  . 
2  3  i  . 
2  34. 

235. 

236. 
2  37. 

238. 

239. 
2  40. 

2  41  . 

2  42  . 
243. 

2  44. 
245. 


DOUBLE  PRECISION  FUNCTION  RM (K  ,GS ,NACT , I , R ,TVAL , IDF1 . PB1 ) 

INTEGER  GS  , GSJ  , GSJL 

REAL*8  C ,  PJL 

KNACT-K-NACT 

ll-I-l 

NACTI-NACT-I 
KNACT1 “KNACT-1 
D«SQRT(PB1 )*R 
L1-MIN0(GS  ,1+1) 

KM -0.  DO 
DO  10  J-l  ,L1 
Jl-J-1 
GSJ-GS-J 

L2 “MI NO(GSJ ,NACTI)+1 
DO  2  0  L-l  ,L2 
GSJL-GSJ-L+1 

PJL-C(I  ,J1 )*C(NACTI  ,L-1 ) *C (KNACT1  ,GSJL) 

DP-D*(L-J) 

CALL  MDTN(TVAL , IDF1 ,DP,A1 ,IER) 

CALL  MDTN(TVAL ,IDF1  ,-DP  ,A2  ,IER) 

2  0  RM-RM  +  PJL*(2  .D0-A1-A2) 

10  CONTINUE 

RM-RM/C(K-1  ,GS-1  ) 

RETURN 

END 

REAL  FUNCTION  D1 (GS  ,G , PG  ,R ,CV) 

INTEGER  GS,G  ,G1  ,S1 
REAL  CV( 4 ) 

Gl-G-1 
D1  -0. 

DO  1  0  J«1  , G 

CALL  MDBIN( J-l  ,G1  ,PG  ,PS,PJJ  ,IER) 

Sl«GS*J+l 

N2  «S1 +4-M0D ( S 1  ,4) 

IDF2-N2-S1 

DP-SQRT(FL0AT(N2  ))*R 

CALL  MDTN(CV(IDF2  )  ,IDF2  ,DP,Al  ,IER) 

CALL  MDTN(CV( I DF2  )  ,IDF2  ,-DP,A2  ,IER) 

10  D1-D1+PJJ*(2  .-A1-A2  ) 

RETURN 

END 


tion  pages  were  taken  from  the 


LIB-0008,  June  1980. 


IMSL  ROUTINE  NAME  -  MOB IN 


-  BINOMIAL  PROBABILITY  DISTRIBUTION  FUNCTION 

-  CALL  MDBIN  (K ,N , P , PS , PK , IER) 

K  -  NUMBER  OF  SUCCESSES  (INPUT) 

N  -  NUMBER  OF  BERNOULLI  TRIALS  (INPUT) 

P  -  PROBABILITY  OF  SUCCESS  ON  EACH  TRIAL  (INPUT) 

PS  -  PROBABILITY  THAT  THE  NUMBER  OF  SUCCESSES  IS 

K  OR  LESS  (OUTPUT) 

PK  -  PROBABILITY  THAT  THE  NUMBER  OF  SUCCESSES  IS 
EXACTLY  K  (OUTPUT) 

IER  -  ERROR  PARAMETER.  (OUTPUT) 

TERMINAL  ERROR 

IER  »  129  INDICATES  THAT  K  WAS  SPECIFIED 
LESS  THAN  ZERO  OR  GREATER  THAN  N. 

IER  -  130  INDICATES  THAT  P  WAS  SPECIFIED 
GREATER  THAN  1  OR  LESS  THAN  0. 

PRECISION/HARDWARE  -  SINGLE/ALL 

REQD.  IMSL  ROUTINES  -  UERTST , UGETIO 

NOTATION  -  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 

Algorithm 

MDBIN  computes  the  probability,  PK,  of  obtaining  exactly  K  successes 
and  PS,  the  probability  of  obtaining  K  or  less  successes  in  N  Bernoulli 
trials,  where  P  is  the  probability  of  success  (Q=l-P  is  probability 
of  failure)  on  each  trial. 

The  following  recursive  relationship  is  used 
Pr(X»J)»Pr(X»J-l)  •  (N+l-J)  ’P/U-Q) 


PURPOSE 

USAGE 

ARGUMENTS 


Then 


Pr(X  not  greater  than  K) 


K 

Z  Pr(X=J) 
J-0 


To  eliminate  the  possibility  of  underflow,  PS  and  PK  are  computed  for- 
ward  from  0  if  K  is  not  greater  than  N*P  and  backward  from  N  otherwise. 
The  scaled  quantities  PS  and  PK,  are  built  recursively  using  an  initial 
value  of  EPS,  where  EPS  is  the  smallest  positive  machine  number.  These 

N 

quantities  are  rescaled  by  a  multiplicative  factor  of  Q  *EPS  if  forward 

N 

computation  is  done  and  P  *EPS  if  backward  computation  is  performed. 


For  the  special  case  P«0 


PK 


if  Kj*0 
if  K-0 


•  9 


IMSL  ROUTINE  NAME  -  MDSTI 

PURPOSE  -  INVERSE  OF  A  MODIFICATION  OF  STUDENTS  T 

PROBABILITY  DISTRIBUTION  FUNCTION 

USAGE  -  CALL  MDSTI  (Q,F,X,IER) 

ARGUMENTS  Q  -  INPUT  PROBABILITY  IN  THE  EXCLUSIVE  RANGE 

(0,1).  (THE  SUM  OF  THE  AREAS  (EQUAL)  IN 
BOTH  TAILS  OF  THE  T  DISTRIBUTION.) 

F  -  INPUT  DEGREES  OF  FREEDOM  FOR  T  DISTRIBUTION. 

F  MUST  BE  GREATER  THAN  ZERO. 

X  -  OUTPUT  VALUE  SUCH  THAT  THE  PROBABILITY  OF  THE 

ABSOLUTE  VALUE  OF  T  BEING  GREATER  THAN  X  IS 

Q. 

IER  -  ERROR  PARAMETER.  (OUTPUT) 

WARNING  ERROR 

IER  -35  INDICATES  OVERFLOW  WOULD  HAVE 
OCCURRED.  X  IS  SET  TO  MACHINE  INFINITY. 
TERMINAL  ERROR 

IER  ■  129  INDICATES  THAT  F  (DEGREES  OF 
FREEDOM)  IS  LESS  THAN  OR  EQUAL  TO  ZERO. 

IER  -  130  INDICATES  THAT  Q  IS  OUT  OF  THE 
EXCLUS IVE  RANGE  (0,1). 

IER  *  131  INDICATES  THAT  AN  ERROR  OCCURRED 
IN  IMSL  ROUTINE  MDNRIS,  THE  INVERSE  NORMAL 
PROBABILITY  DISTRIBUTION  FUNCTION. 

IER  ■  132  INDICATES  THAT  CONVERGENCE  WAS 
NOT  ACHIEVED.  THIS  CAN  ONLY  OCCUR  FOR 
F  LESS  THAN  2.0  . 

PRECISION/HARDWARE  -  SINGLE/ALL 

REQD.  IMSL  ROUTINES  -  MDNRIS , MERFI , UERTST , UGETIO 

NOTATION  -  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 


REMARKS 


NOTE  THAT  MDSTI  DOES  NOT  PROVIDE  THE  ACTUAL  T  INVERSE. 
FOR  P  EQUAL  TO  THE  PROBABILITY  THAT  A  STUDENTS  T 
RANDOM  VARIABLE  IS  LESS  THAN  X,  THAT  INVERSE  CAN  BE 
OBTAINED  BY  THE  FOLLOWING  RULES. 

A.  FOR  P  IN  THE  RANGE  (0.0, 0.5),  CALL  MDSTI  WITH 
Q  ■  2*P  AND  NEGATE  THE  RESULT  X. 

B.  FOR  P  IN  THE  RANGE  (0.5, 1.0),  CALL  MDSTI  WITH 
Q  -  2M1-P)  . 


Algorithm 


This  subroutins  computes  X  such  that 

X 

Q  -  1-  /  fi(t)dt 

-X  F+l 

-204 


here  f  (t) 


1  +  r 


for  -ao  <  t  < 


See  the  following  figure. 


See  reference: 

1.  —Hill,  G.W.,  "Algorithm  396,  Student's  t-quantiles" ,  CACM,  13(10)1970, 
619-620. 


2.  Roe,  G.M. ,  Programs  for  the  Incomplete  Beta  and  Gamma  Functions 
and  their  Inverses ,  General  Electric  Technical  Information 
Series ,  January,  1969. 

Accuracy 

Test  results  indicate  5  significant  digit  accuracy.  In  addition.  Table 
26.10,  p.  999  of  the  Handbook  of  Mathematical  Functions ,  edited  by 
M.  Abramowitz  and  I. A.  Stegun,  was  duplicated  by  MDSTI. 

Example 

In  this  example  X  is  calculated  such  that  the  probability  of  the  ab¬ 
solute  value  of  a  random  variable  following  Student's  t  distribution 
with  two  degrees  of  freedom  being  greater  than  X  is  0.8. 

Input : 

INTEGER  IER 
REAL  Q,F,X 
Q  -  0.8 

F  «  2. 

CALL  MDSTI  (Q,F,X,IER) 

Output : 


X 


0.28867 


IMSL  ROUTINE  NAME 


MDTN 


PURPOSE 


NON-CENTRAL  T  PROBABILITY  DISTRIBUTION 
FUNCTION 


USAGE 


CALL  MDTN  (TVAL, IDF,D,P,IER) 


ARGUMENTS  TVAL 
IDF 
D 
P 


IER 


INPUT  VALUE  TO  WHICH  INTEGRATION  IS  PERFORMED. 
INPUT.  INTEGER  DEGREES  OF  FREEDOM. 

INPUT.  NON-CENTRALITY  PARAMETER. 

OUTPUT.  PROBABILITY  THAT  A  RANDOM  VARIABLE 
DISTRIBUTED  AS  T  WITH  NON-CENTRALITY 
PARAMETER  D  IS  LESS  THAN  OR  EQUAL  TO  TVAL. 
ERROR  PARAMETER.  (OUTPUT) 

TERMINAL  ERROR 

IER  -  129,  INDICATES  DEGREES  OF  FREEDOM  OUT 
OF  RANGE.  IDF  MUST  BE  GREATER  THAN  ZERO. 


PRECISION/HARDWARE  -  SINGLE /ALL 


REQD.  IMSL  ROUTINES  -  MDNOR, MDTNF , MERRC-ERFC , UERTST , UGETIO 

NOTATION  -  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 


Algorithm 

MDTN  calculates  the  probability  P  that  a  random  variable  Z  which  follows 
the  non-central  t  distribution,  with  non-centrality  parameter  D,  is  less 
than  or  equal  to  TVAL. 

Letting  f(t)  be  the  non-central  T  probability  density  function  with 
non-centrality  parameter  D,  the  indicated  P  value  is  returned. 


See  reference: 

Owen,  D.B.,  "A  special  case  of  the  bivariate  non-central  t-distribution" 
Biometrika,  52,  1965,  437-446. 


Example 

This  example  illustrates  use  of  MDTN  with  a  positive  non-centrality 
parameter,  D*6.0. 
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